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ABSTRACT 

The cosmic dark ages ended a few hundred milhon years after the Big Bang, when the first stars 
began to fill the universe with new light. It has generally been argued that these stars formed in 
isolation and were extremely massive - perhaps 100 times as massive as the Sun. In a recent study, 
Clark and collaborators showed that this picture requires revision. They demonstrated that the 
accretion disks that build up around Population HI stars are strongly susceptible to fragmentation 
and that the first stars should therefore form in clusters rather than in isolation. We here use a 
series of high-resolution hydrodynamical simulations performed with the moving mesh code AREPO 
to follow up on this proposal and to study the influence of environmental parameters on the level 
of fragmentation. We model the collapse of five independent minihalos from cosmological initial 
conditions, through the runaway condensation of their central gas clouds, to the formation of the first 
protostar, and beyond for a further 1000 years. During this latter accretion phase, we represent the 
optically thick regions of protostars by sink particles. Gas accumulates rapidly in the circumstellar 
disk around the first protostar, fragmenting vigorously to produce a small group of protostars. After 
an initial burst, gravitational instability recurs periodically, forming additional protostars with masses 
ranging from ~ 0.1 to 10 M©. Although the shape, multiplicity, and normalization of the protostellar 
mass function depend on the details of the sink-particle algorithm, fragmentation into protostars with 
diverse masses occurs in all cases, confirming earlier reports of Population HI stars forming in clusters. 
Depending on the efficiency of later accretion and merging. Population HI stars may enter the main 
sequence in clusters and with much more diverse masses than are commonly assumed. 
Subject headings: cosmology: theory — early universe — hydrodynamics — methods: numerical — 
stars: formation 



1. INTRODUCTION 

Over the past decade, a consensus has emerged on how 
the first, so-called Population HI (Pop HI), stars may 
have formed o ut of the hydrogen and helium forged in 
the Big Bang ([Ba rkana fc Loebl [20011 : iBromm fc LarsonI 
120041 : iBromm et a l. 2009). Non-linear structure forma- 
tion from an initially nearly featureless universe began 
when dark matter minihalos with masses of order 10^ M© 
collapsed at redshifts z ^ 20 — 50, confining gas within 
their gravitational potential wells (jHaiman et al.l [19961 : 
[Tegmark et al. 19^3). Cooling through ro- vibrational 
transitions of newly formed molecular hydrogen then 
triggered runaway collapse at the center of the halo, 
resulting in the formation of a protostar with density 
more than twe nty orders of magnit ude higher than the 
cosmic mean (jYoshida et al.l |2008[ ). The envelope of 
the accreting protostar remained hot due to the lack 
of cooling by metals and dust, so that accretion rates 
were on average higher than in star-formation regions to- 
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dav (iMac Low fc Klessedl2004l : iMcKee fc Ostriker[l2007l : 
IZinnecker fc Yorke 2007^ With a few exceptions 
(|Turk et~l.. 2009), simulations of thi s initial collapse 
phase have shown no fragmentation (lAbel et al.l [2QQ2[: 
Bromm et al[l2002l: [B omm fc Loebll2004l:lYoshida et al. l 



20061 : 10^ Shea fc NormanI l2007l : lYoshida et al. I I2008D . 

leading to the conclusion that the first stars formed in 
isolation and were extremely massive. 

In contrast, studies of present-day star formation have 
generally found fragmentation to occur shortly after the 
formation of the f ir st protostar ([Klessen fc Burkerlj[2000l: 
Bate et all 120031 : iKrumholz et al.l 120091 : [Peters et aP 



20T3). Following up on this result, sink particles were 



recently used in studies of the fragmentation of primor- 
dial gas in minihalo s (Clark et al. 2008; Stacy et al. 20Tq[ : 
[Clark et al.ll2011a[ ). They found that the metal- free gas 
clouds fragment strongly, with the details of the process 
depending on the degree of turbulence in the halo. Fo- 
cusing on the dynamical evolution of the high-density gas 
in the central regions of a minihalo, Clark et al. (201 13) 
demonstrated that the protostellar disks around primor- 
dial stars accrete from the infalling envelope faster than 
they can transfer their mass onto the central object. As a 
result, they rapidly become gravitationally unstable and 
fragment to build-up binary or higher-order multiple stel- 
lar systems. These results challenge the idea that the first 
stars formed in isolation and give rise to a number of new 
questions: What is the mass spectrum of Pop III stars 
in groups? Is it different from the previously preferred 
single mode of star formation? How does it depend on 
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TABLE 1 
Simulation parameters 



Simulation 


Size [kpc] 


Particles 


Mdm [M©] 


Mdm,ref [M©] 


Mgas [Mo] 


cr8 


Mvir [Mo] 


r^riv [pc] 


A 


^coll 


MH-1 


1000 


512^ 


272 


3.53 


0.72 


0.81 


5.8 X 10^ 


150 


0.059 


18.6 


MH-2 


500 


256^ 


272 


3.53 


0.72 


0.9 


3.0 X 10^ 


110 


0.055 


19.5 


MH-3 


250 


128^ 


272 


3.53 


0.72 


1.2 


2.3 X 10^ 


94 


0.073 


20.9 


MH-4 


500 


256^ 


272 


3.53 


0.72 


1.1 


3.1 X 10^ 


97 


0.044 


22.6 


MH-5 


500 


256^ 


272 


3.53 


0.72 


1.3 


1.8 X 10^ 


58 


0.038 


31.7 



The comoving box sizes, particle numbers, initial DM masses, refined DM masses, gas masses, and normalizations used in the 
simulations, as well as the viral masses, virial radii, spin parameters, and collapse redshifts of the first minihalos t hat form. 
The halo p roperties agree well with the results of previous studies (e.g., iMachacek et al.h2001UYoshida et al.ll2003i : IGao et all 
[20071 : fO^ ea &; Norman 2007). 



the dynamical characteristics of the host halo? How does 
it influence subsequent cosmic evolution? 

To address these issues, we here report on a se- 
ries of hydrodynamical simulations performed with the 
quasi-Lagrangian moving mesh code AREPO (jSpringell 
[20103) . Its high accuracy, flexibility, and efficiency 
allow us to follow the evolution of the gas in five 
statistically independent minihalos from fully cosmo- 
logical initial conditions down to the optically thick 
regime where primordial protostars are born. We use 
a sink-particle approach that is similar to previously 
employed techni ques to capture the subsequent accre- 
tion phase (e .g., Bate et al.Ml995l; iKrumho lz erall 120041 : 
iJappsen et al . 2005; Federrath et al. 2010), which allows 
us to reach well beyond the formation of the first pro- 
tost ar, and circumvent the limitations posed in previous 
high-resolu tion ab initio calcula tions of primordial star 
formation ([Yoshida et al.l [20081 ). The five realizations 
give us an indication of how common fragmentation is, 
and of the shape of the mass function of Pop III proto- 
stars during the early stages of accretion. 

The structure of our work is as follows: In Section 2, 
we describe the numerical setup and physical ingredients 
of the simulations. In Section 3, we present the results 
of the simulations, followed by a discussion of radiation 
feedback, a resolution study, and the caveats of the sink 
particle treatment. Finally, in Section 4 we summarize 
our results and draw conclusions. All distances quoted 
in this paper are in proper units, unless noted otherwise. 

2. NUMERICAL METHODOLOGY 

We here provide details about the set-up of the pure 
dark matter (DM) simulations, the employed resimula- 
tion technique, and the main simulation runs with the hy- 
drodynamic moving mesh code AREPO. We also discuss 
the implementation of the on-the-fly mesh refinement 
technique, the sink-particle treatment, and the chemistry 
and cooling model. 

2.1. Dark Matter Simulations 

The DM-only simulations are initialized in cosmolog- 
ical boxes with 250, 500 and 1000 kpc (comoving) on a 
side, and employ 128^, 256^ and 512^ particles of mass 
^dm — 272 M©, so that minihalos with a characteristic 
mass of 5 X 10^ Mq are resolved by 2 x 10^ particles and 
a comoving gravitational softening length of 68 pc. We 
use an initial fluctuation power spectrum for a A cold 
dark matter (ACDM) cosmology at a starting redshift 
of z = 99 with matter density = 1 — = 0.27, 
baryon density Qi, = 0.046, Hubble parameter h = 



i^o/100 kms~^ Mpc~^ = 0.71 (where Hq is the present 
Hubble expansion rate) and a spectral index = 0.96. 
These paramet ers are based on the five-year WMAP re- 
sults (Komatsu et al.l [2009). We use a range of values 
for the linearly extrapolated present-day normalization 
as (as given by Table 1), which allows us to mimic rare 
overdense regions of the universe that cannot be captured 
with the limited box sizes employed here. 

All our DM simulations were evolved with the 
TreePM/SPH code GADGET-3 (last described in 
■Springel 2005) until the first minihalo with a virial mass 
exceeding 5 x 10^ Mq collapses, thereby providing a tar- 
get region for our subsequent resimulations. In this 
study, we restrict ourselves to the analysis of five inde- 
pendent realizations of the minihalos identified in the 
large-scale DM runs, denoted MH-1 to MH-5. The most 
important parameters of the simulations and the prop- 
erties of the minihalos are summarized in Table 1. 

2.2. Resimulation Technique 

Once the location of a minihalo in the parent simula- 
tion has been determined, we construct a high-resolution 
study of this object by first identifying all particles within 
the virial radius and a sufficiently large boundary region 
around it. These particles are then traced back to their 
initial conditions, yielding the Lagrangian region that 
formed the object. To construct new initial conditions, 
each particle in this region is replaced by 64 DM par- 
ticles and 64 mesh-generating points which produce the 
space-filling cells that are used for our hydrodynamic cal- 
culations. 

Cells and DM particles outside of the high-resolution 
region are replaced by ever higher mass particles with 
increasing distance from the target halo, such that the 
resolution coarsens significantly far away from the tar- 
get halo, yet the gravitational tidal field that influences 
its formation is still followed with high accuracy. This 
step reduces the total number of particles and cells in the 
refined initial conditions to 2 x 10^. The simulation 
box is then centered on the target halo and reinitialized 
at z = 99 based on the original realization of the density 
field, but augmented with additional small-scale power in 
the high-resolution region that can now be represented. 
The initial DM particle and cell masses in the high- 
resolution region before any further run-time refinement 
are given by Mdm,ref = (1 - ^^6/^m)Mdm/64 3.53 M© 
and Mgas = (^^6/^^m)^dm/64 0.72 M0, respectively. 
We use a comoving gravitational softening length of 17 pc 
for the refined DM component. 
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Fig. 1. — A test of the classical 'Truelove et al.l (|1998l ) criterion in AREPO for the lBoss &: Bodenheimerl (|1979l ) isothermal collapse problem. 
The individual panels compare simulations with on-the-fly mesh refinement with 1, 2, 4, and 8 cells per Jeans length, respectively. We 
show the density-squared weighted mass density projected along the line of sight in a box with 200 AU on a side, centered on one of the 
two main clumps. Similar to the result found in AMR codes, the amount of artificial fragmentation increases dramatically once the local 
Jeans length is resolved by less than four cells. 



2.3. The Moving Mesh Code AREPO 

We follow the collapse of the gas in the refined mini- 
halos with the c osmological moving mesh code AREPO 
(ISpringell [20103) . AREPO is a second-order accurate 
finite volume method that solves the Euler equations 
based on a piece-wise linear reconstruction and the cal- 
culation of hydrodynamical fiuxes at every cell face with 
an exact Riemann solver. The principle difference of 
AREPO compared to Eulerian mesh codes is that its 
computational mesh is constructed as the Voronoi tessel- 
lation of a set of mesh-generating points. These points 
can be moved with the fiow velocity itself, making the 
mesh automatically adaptive in a Lagrangian fashion. 
This technique greatly reduces the numerical diffusiv- 
ity of mesh-based hydrodynamics, especially when large 
bulk fiows are present. In fact, the results of the code 
become fully Galilean-invariant, whereas the truncation 
error of ordinary Eulerian mesh codes depends on the 
bulk velocity of the system. In addition, the unstruc- 
tured Voronoi mesh of AREPO avoids the introduction 
of preferred directions, which are present in Cartesian 



meshes. 

The novel AREPO scheme hence combines the accuracy 
of mesh-based hydrodynamics with the natural adaptiv- 
ity and translational invariance usually only provided by 
the smoothed particle hydrodynamics (SPH) technique. 
In terms of hydrodynamical accuracy, the grid-based ap- 
proach of AREPO alleviat es a number o f shortcomings 
encountered with SPH ( Mo naglianl l2005f ) . Among these 
are the inherent noise of the kernel estimates, the ar- 
tificial viscosity, and the slow convergence rate of SPH 
in three dimensions (|Springell l2010b[ ) . Further impor- 
tant improvements of AREPO lie in the more accurate 
treatment of shocks a nd turbulence, and of fiuid in- 
stabilities (Agertz et al. 2007). Compared to adaptive 
mesh refinement (AMR) codes, a significant advantage of 
AREPO lies in its ability to continuously adjust its reso- 
lution when density fiuctuations grow under self-gravity. 
This key feature of Lagrangian codes is ideal for gravita- 
tional collapse problems. If the bulk velocities are large, 
AREPO can use larger timesteps than AMR codes and 
exhibits lower numerical diffusion errors at comparable 
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Fig. 2. — A sequential zoom-in on the gas in each of the five minihalos. The panels show the density-squared weighted number density of 
hydrogen nuclei projected along the line of sight. The gas virializes on a scale of ~ 5 kpc (comoving), followed by the runaway collapse of 
the central ~ 1 pc, where the gas becomes self-gravitating and decouples from the dark matter. In the final stages of the collapse, a fully 
molecular core on a scale of a few hundred AU forms. The visible turbulence induced by the virialization of the dark matter halo survives 
down to the smallest scales and later influences the fragmentation of the gas. 



spatial resolution. 

We note that the computational speed of AREPO in 
simulations of cosmic structure formation is roughly on 
par with GADGET-3 for an equal number of hydro- 
dynamic resolution elements, despite its comparatively 
complex mesh-construction calculations. This is in part 
due to the large cost of accurate calculations of self- 
gravity, which dominate over the hydrodynamics both 
in AREPO and GADGET-3. Both codes compute the 
gravitational field with a fully adaptive combination of a 
particle-mesh (PM) solver for large-scale forces, and a hi- 
erarchical multipole expansion for short-range forces. In 
any case, because of the higher accuracy of AREPO com- 
pared to SPH for an equal number of resolution elements, 
this means that AREPO is the more efficient method. 



The faster convergence rate of AREPO also implies that 
its accuracy advantage increases further if more resolu- 
tion elements are employed. 

The mesh-generating points of AREPO can still be in- 
terpreted as Lagrangian fiuid particles, if desired, which 
greatly simplifies the (re) use of well-tested techniques 
from SPH codes, such as implementations of gas chem- 
istry or sink particles, an approach we employ extensively 
in this study. Also, even though AREPO 's fiuid cells fol- 
low the motion of the gas and hence maintain a roughly 
constant gas mass, it is sometimes necessary to refine the 
mass resolution per cell during a calculation. This can be 
done accurately in the mesh-based approach of AREPO, 
as we discuss next. 
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Fig. 3. — See Figure 2 for caption. 
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2.4. Mesh Refinement 

It is well established that an essential prerequisite for 
reliable hydro dynamic collapse simulations lies in a nu- 
merical resolution of the Jeans length. In SPH simula- 
tions, the Jeans mass associated with each particle must 
exceed the mas s within the smo othing length by a factor 
of about two ([Bate fc Burkerti il997i) , while grid-based 
codes rely on the so-ca lled ^Truelove c riterion' for satis- 
fying this requirement (|Truelove et al.| [l998). which stip- 
ulates that in an isothermal collapse calculation at least 
four cells per Jeans length should be used. Collapse sim- 
ulations may therefore require additional run-time refine- 
ment to ensure that the Jeans length is resolved and ar- 
tificial fragmentation is avoided. 

The addition of cells in grid codes is comparatively 
straightforward, while special care must be taken in SPH 
simulations t o avoid substantial noise when new particles 
are inserted (|Kitsionas fc WhitwortE 12002). In AREPO, 
a very simple and robust approach exists: if a cell ful- 
fils a predefined refinement criterion, it is split by in- 
serting a further mesh-generating point at the position 
of the generator of the original cell, displaced by a ran- 
dom offset that is very small compared to the size of the 
cell. The mass, momentum and energy of the original 
cell are then distributed conservatively among the two 
new cells, weighted by their respective volumes, keeping 
the density, velocity and pressure of the original cell un- 
changed. Over the course of a few timesteps, the code 
then separates the mesh-generating points of the two 
cells and moves them closer to their geometric centers 
of mass, as a result of the mesh-steering corrections em- 
ployed by AREPO to maintain a regular structure of the 



mesh. We note that the numerical fiuxes required for this 
mesh readjustment arise consistently from the Riemann 
solutions calculated by the code, as these take the mesh 
motion fully into account. 

We here briefly investigate how the resolution of 
the Jeans length affects the amount of artificial 
fragmentation found in AREPO, using a series of 
isothermal collapse sim u lation s in the formulation of 
iBurkert fc Bodenheimerl (Il993f), which is based on a 
test initially introduced bv lBoss fc Bodenheimerl (|1979[ ). 
This setup consists of a 1 M© cloud with radius r = 
5 X 10^^ cm and sound speed Cg = 1.66 x lO^cms"^ 
in solid-body rotation with an angular velocity oi uj = 
7.2 X 10~^^rads~^ and an m = 2 density perturbation 
of the form 



=po[l + 0.1cos(2(/))] , 



(1) 



where po = 3.82 x 10~^^gcm~^ is the underlying con- 
stant density and (j) is the azimuthal angle around the 
rotation axis. We refine the gas whenever the local Jeans 
number, here defined as the radius of a cell divided by 
the local Jeans length, increases above a predefined value. 
We estimate the cell radius as = (3y/47r)^/^, where V 
is the volume of the cell. Since the mesh-steering mo- 
tions ensure that the cells do not become too distorted, 
this provides a good estimate of the size of a cell. In Fig- 
ure 1, we show the state of the gas in one of the two main 
clumps for a minimum of one, two, four, and eight cells 
per Jeans length. For isothermal gas, artificial fragmen- 
tation sets in above a Jeans number of approximately 
1/4, which is similar to the result found in AMR codes 
(jTruelove et al.|[l998[ ). However, this should be consid- 
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Fig. 4. — Temperature versus number density of hydrogen nuclei in the fully cosmological simulations. The mass in each bin is color-coded 
from blue (minimum) to yellow (maximum). The evolution of the gas is similar to the results of previous studies, with the exception that 
HD cooling becomes important in simulations MH-2 and MH-3 and leads to a prolonged cooling phase at nu > 10^ cm~^. 



ered a lower limit on the required resolution. 

In the present study, we activate the refinement crite- 
rion in our cosmological simulations above a density of 
nn = 1 cm~^ in the central 200 pc of the halo, which 
is larger than the maximum virial radius of the mini- 
halos investigated here (see Table 1). We ensure that 
the local Jeans length is resolved by at least 128 cells 
until the number density of hydrogen atoms exceeds 
nH = lO^cm"^, at which point we extract the cen- 
tral 1 pc of the simulations and deactivate further refine- 
ment. We reinitialize the simulations in this smaller box 



with refiective boundary conditions (see lSpringelll2010al ). 
which ensures that the gas remains confined by exter- 
nal pressure. We also discard all DM particles, since the 
gravitational potential at these densities is dominated 
by the gas: the radius at which the gas mass exceeds 
the DM mass by an order of magnitude is 0.76, 0.67, 
0.58, 0.71 and 1.1 pc, respectively, which is comparable 
to the spatial extent of the resimulations. Due to the de- 
activation of the refinement, the particle mass within the 
central 1000 AU remains roughly constant at 10~^ M©. 
We note that shocks from infiows into the innermost re- 
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Fig. 5. — Same as Figure 4, but for the follow-up high resolution simulations. In simulations MH-2 and MH-3, the gas heats up more 
violently after the prolonged cooling phase due to HD cooling, which creates shocks that are visible as 'fingers' of hot, underdense gas. At 
densities greater than ~ lO"*^^ cm""^, the gas becomes optically thick to cooling radiation and shortly thereafter forms a protostar. 



gion are still well resolved under these conditions, thanks 
to the ability of AREPO to capture shocks on ~ 1 — 2 
cells. Assuming a temperature of 1000 K, the Truelove 
criterion is thus expected to be violated at a density 
of lO^^cm"^, which is approximately the density at 
which the sink particles are inserted. 

2.5. Sink Particles 

The evolution of the gas beyond the runaway col- 
lapse of the first protostellar core is captured with a 
sink-particle algorithm that exploits the hybrid nature 



of AREPO. Since the mesh-generating points are La- 
grangian tracers of the fluid, sink particles may be 
treated very similarly t o implementations commonly 
used i n SPH simulations (iBate et al.lll995l : iBromm et al.l 
I2OO2I : iJappsen et al.l l2005f K However, instead of using 
a density threshold, we identify candidate sink particles 
whenever the Jeans number of a cell falls below a pre- 
defined critical Jeans number JcHt = 1/8- For a sink 
particle to form, the gas within rgink = h/ Jcvit is further 
required to be bound and have a negative velocity diver- 
gence. All mesh-generating points within rgink are then 
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Fig. 6. — The panels on the left show the spherically averaged number density of hydrogen nuclei and enclosed gas mass as a function of 
radius just before the formation of the first protostar. The panels on the right show the spherically averaged specific angular momentum, 
temperature, and velocity dispersion in units of the sound speed as a function of enclosed gas mass. The density, enclosed gas mass, and 
angular momentum profiles are very similar overall, while the thermal evolution of the gas displays some scatter due to the activation of HD 
cooling in simulations MH-2 and MH-3. In these two minihalos, the gas heats up later but more violently after becoming gravitationally 
unstable. The velocity dispersion shows no convincing correlation with the thermal history of the gas and is always close to Mach numbers 
M ~ 1, indicating transsonic turbulence. An interesting feature is the simultaneous collapse of a second gas cloud in simula tion MH-2, 
which is visible as a density peak at about 1000 AU from the primary cloud. Similar behavior was found in another recent study ([Turk et al.l 
[20091 ). 



removed and replaced by a collisionless sink particle with 
a gravitational softening length of racc/3, where race is 
a predefined accretion radius that is set independently 
from the initial sink particle radius. In practice, only 
200 mesh-generating points are removed during this 
step, since the spatial resolution around the sink particle 
decreases at larger radii. The sink particle is placed at 
the center of mass of the removed cells with a velocity 
determined by linear momentum conservation, while the 
angular momentum and internal energy of the gas are 
discarded. The additional factor 1/3 in the gravitational 



softening is used to avoid artificial fragmentation, which 
might occur if gravitational forces on the gas are reduced 
on the scale of the accretion radius. 

Accretion onto existing sink particles occurs if the 
mesh-generating point associated with a candidate cell 
falls within the accretion radius race of the sink particle 
it is most bound to. This method exploits the Lagrangian 
nature of the mesh-generating points, and yields the (in- 
cremental) rate at which mass flows onto sink particles. 
We have also tested an implementation with a more strin- 
gent criterion, where in addition to being bound to the 
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Fig. 7. — The formation of a protostellar cluster at the center of the minihalos using standard sink particles. The panels show the density- 
squared weighted number density of hydrogen nuclei projected along the line of sight. Black dots and crosses denote protostars with masses 
below and above 1 M©, respectively. The process of initial disk formation and fragmentation is remarkably similar in all minihalos (see also 
I Clark et al . 2011b), after which N-body effects become important and lead to relatively unique configurations. For example, in simulation 
MH-4 dynamical interactions have led to the ejection of a low-mass protostar after only ~ 50 yr. This occurs significantly later in the other 
four minihalos. 
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Fig. 8. — The central 2000 AU after 1000 yr of continued fragmentation and accretion. Black dots, crosses and stars denote protostars with 
masses below IM0, between 1 and 3M0, and above 3M0. A relatively rich protostellar cluster with a range of masses has survived 
in each case. In a few minihalos, low-mass protostars have been ejected out of the central gas cloud, such that they are no longer visible 
here. In simulation MH-2, two independent clumps have collapsed almost simultaneously and formed their own clusters before eventually 
merging (see also Figure 6). 
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Fig. 9. — The mass accretion histories of the protostars (left panels), and the average accretion rates of the entire ensemble of protostars 
(right panels) using standard sink particles. In the left panels, each line denotes the evolution of an individual protostar. The thick dashed 
lines in both panels denote the cumulative values. Mergers are indicated by instantaneous jumps in mass. We find that star formation is 
not restricted to a single burst, but occurs continually for the entire simulated timespan. In every minihalo, between 5 and 15 protostars 
with masses ranging from 0.1 to nearly 10 are formed. The total accretion rates are nearly constant over time at a few 0.01 Mq yr~^. 



sink particle, the semimajor axis of the respective two- 
body system must fall below the accretion radius. How- 
ever, we have found that these additional criteria made 
no appreciable difference in practice. 

Once a cell has been accreted, its mass and momen- 
tum is added to the sink particle and the correspond- 
ing mesh-generating point is removed. After this step 
has been performed for all candidate cells, the mesh 
is reconstructed such that the volume associated with 
the removed cells is distributed among the remaining 
cells around the sink particle. As conserved quantities 
are weighted with the new volumes, this tends to make 



their densities and pressures artificially small. In the 
last section, we present a resolution study with varying 
accretion radii to show that this caveat artificially re- 
duces the amount of fragmentation, and increases the 
typical fragment mass. As our fiducial accretion radius, 
we choose a value of race = 100 R©, which is close to 
the maximum physical size of accreting Pop III stars 
(Hosokawa & Omukai" 2009). 

Mergers between sink particles occur whenever the to- 
tal energy of the respective two-body system is negative 
and the semimajor axis falls below the accretion radius. 
This is motivated by our above choice, which yields an 
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Fig. 10. — Same as Figure 9, but using adhesive sink particles. In this case fewer protostars with systematically higher masses are formed, 
although the total amount of gas within protostars is nearly identical. 



approximate upper limit on the physical size of the pro- 
tostars. In cases where more than one sink particle lies 
within the accretion radius of another sink particle, the 
pair with the highest binding energy merges. For com- 
parison, we have also implemented an alternative merg- 
ing criterion, where sink particles merge whenever their 
separation drops below race, and no further checks on 
their energies are enforced. The purpose of this formula- 
tion is to maximize the efficiency of merging between pro- 
tostars, since we do not capture the gasdynamical friction 
between real protostars. 



2.6. Chemistry and Cooling 



The chemical and thermal evolution of the gas is mod- 
eled as described in IClark et all (|2Qlla[ ) . This network 
consists of 45 reactions amongst twelve chemical species: 
H, H+, H-, H+, H2, He, He+, He++, D, D+, HD, and 
free electrons. We account for all relevant cooling pro- 
cesses, including electronic excitation of H, He and He+, 
cooling from the recombination of H+ and He+ , Compton 
cooling, and bremsstrahlung. In practice only a few pro- 
cesses play an important role at the densities and temper- 
atures investigated here - H2 rotational and vibrational 
line cooling, H2 collision-induced emission (CIE) cool- 
ing, H2 collisional dissociation coolin g, and heatings; d ue 
to three-body H2 formation (see also lTurk et al.|[2QlTh . 

Collisional dissociation of H2 is modeled with the rate 
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Fig. 11. — The protostellar mass function after 1000 years of continued fragmentation and accretion. The dark and hght shadings 
distinguish the mass functions obtained for standard and adhesive sink particles, respectively. Despite very aggressive merging, a small 
cluster of protostars with a range of masses is formed even in the latter case. In the bottom right panel, we also show the cumulative mass 
functions obtained by summing up the contributions from the individual minihalos, and renormalized for better visibility. The resulting 
distribution is relatively flat between ~ 0.1 and ~ 10 M©, indicating that most of the mass is locked up in high-mass protostars. 



given bv lMartin et al.l ([T 



from which the three-body 
formation rate is derived by applying the principle o f 
detailed balance (Flower & H axrisl |2007|; iGloverl [2QQ8h . 
These rates are intermediate in terms of the rates dis- 
cussed in iTurk et aL (2011), and their uncertainty is re- 
flected in the substantial variation of the thermal and 
morphological evolution of primordial gas clouds. In the 
present study, we do not investigate how the collapse 
of the gas is affected by our choice of the three-body 
rate. H owever, we note that the rate used by Clark et al. 
(j2Qllb [). taken from Abel et al. ( 2002 ). is significantly 
smaller than the rate adopted here, and the fact that 
we see qualitatively similar behavior in both simulations 
gives us confidence that our main results do not depend 
to any great extent on the choice of three-body rate, 
although the quantitative details will likely have some 
dependence on the choice of rate. 

The cooling of the gas by H2 lines in the optically thin 
regime is modeled with the low-density cooling rates for 
collisions between H2 molecules and H and He atoms, 
H2 m olecules, protons and electrons (jGlover Abell 
l2008f ). accounting for the transition to local thermo- 
dynamic equilibrium level populations at gas densities 
riH ^ lO^cm"^. At densities above nn ~ lO^cm"^, the 
strongest of the H2 ro- vibrational lines become optically 
thick, reducing the effectiveness of H2 line cooling. To 
account for this effect, we use an approach based on the 
Sobolev approximation ( Yoshida et al...2006>) . We write 
the H2 cooling rate as 



Ah2 = ^A^ul^ul/^e: 



(2) 



where 11^ is the number density of hydrogen molecules in 
upper energy level A^^i is the energy difference be- 
tween this upper level and a lower level A^i is the spon- 
taneous radiative transition rate for transitions between 
u and /, and /3esc,ui is the escape probability associated 
with this transition, i.e. the probability that the emit- 
ted photon can escape from the region of interest. We 
fix flu by assuming that the H2 level populations are in 
local thermodynamic equilibrium, and write the escape 
probabilities for the various transitions as 



Pesc,u\ — 



1 - exp(-rui) 



(3) 



where we use the approximation that 

Tui auiLs , (4) 

where aui is the line absorption coe fficient and Lg is the 
Sobolev length (Yoshida et al. 2006). In the classical, 
one-dimensional spherically symmetric case, the Sobolev 
length is given by 



Vth 



\dvr/dr\ ' 



(5) 



where Vth is the thermal velocity, and dv^/dr is the radial 
velocity gradient. In our three-dimensional simu lations, 
we generalize this as (|Neufeld fc Kaufmanlll993l ) 



'^th 

IV -vl 



(6) 



To prevent the H2 cooHng rate from being reduced by an 
unphysically large amount in regions with small velocity 
gradients, we limit Lg to be less than or equal to the local 
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Fig. 12. — The radial velocity and ratio of the radial velocity to the escape velocity of all protostars as a function of distance to the 
center of mass, shown for standard (left panels) and adhesive (right panels) sink particles. The critical ratio of unity, where protostars are 
assumed to escape from the central gas cloud, is denoted by the dashed line. The escape velocity is determined by using the total mass 
enclosed within the current distance of each protostar from the center of mass. Black dots, crosses and stars denote protostars with masses 
below 1 M0, between 1 Mq and 3 M©, and above 3 Mq. In our standard implementation of sink particles, a number of low-mass protostars 
obtain high radial velocities and escape from the central gas cloud. They stop accreting after they receive substantial radial velocities 
during close encounters with other protostars. For adhesive sink particles, this occurs significantly less often. In this case no protostars 
reside within the central 2^ 100 AU, since momentum conservation acts to move merged protostars to larger annuli. 



Jeans length, Lj. We justify this choice by noting that 
there are strong density gradients in the gas on scales 
L Lj^ and so we expect the bulk of the contribution 
to the H2 line absorption to come from material within 
only a few Jeans lengths. 

At very high densities (nn > lO^^cm"^), CIE cool- 
ing from H2 becomes more effective than H2 lin e cooling 
(|Qmukai fc NishilfTOQSl : iRipamonti fc Abelll2QQ4f ) . We ac- 
count for the reduction of the CIE cooling rate due to the 
effects of continuum opacity using the following prescrip- 



tion (|Ripamonti et alJl2QQ2l : iRipamonti fc Abelll2QQ4f ): 

/I — e~^^^^ \ 

AciE,thick = AciE,thin min , 1 , (7) 

\ rciE J 

where 

/ X 2.8 

Finally, we account for the fact that each time an H2 
molecule is collisionally dissociated, 4.48 eV of thermal 
energy is lost by the gas, while every time that a new H2 
molecule is formed by the three-body process, the gas 
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Fig. 13. — A comparison of the fragmentation of the gas in simulation MH-1 with and without infrared radiation emitted by accreting 
protostars. The panels show the density-squared weighted temperature projected along the line of sight. The radiation slightly heats and 
puffs up the disk, such that the fragments form at somewhat larger distances from th e center. T he qualitative nature of radiation feedback 
is therefore very similar to what studies of present-day star formation have shown (K rumholz et al.t 2009: Pet ers et al.i i2010l. However, the 
effect is significantly reduced here due to the high temperature of the gas and its very efficient cooling by molecular hydrogen lines. 



gains 4.48 eV. 

3. RESULTS 

We here discuss the cohapse of the gas in the mini- 
halos up to the formation of the first protostar, and the 
subsequent fragmentation and accretion that leads to the 
build-up of the protostellar cluster. We then proceed to 
investigate the influence of radiation, present a resolu- 
tion study, and discuss the caveats of the sink-particle 
algorithm. 

3.1. Collapse of Gas in Minihalos 

In Figures 2 and 3, we show the collapse of the gas 
in all five minihalos from cosmological to protostellar 
scales. As gas falls into the DM halo, shocks thermal- 
ize the gravitationally generated infall motions and heat 
the gas to the halo virial temperature. At the same time, 
irregular motions of the dark matter seed subsonic tur- 
bulence in the gas, which cascades down to ever smaller 
scales and later influences the fragmentation of the pro- 
tostellar cloud. Once enough hydrogen molecules have 
formed, the gas at the center of the halo begins to cool, 
becomes gravitationally unstable, and d ecouples from the 
dark matter on a scale of a parsec (|Abel et al.l I2QQ2I : 
iBromm et al.l I2QQ2 ). At significantly higher densities, 
three-body reactions render the ^ as fully molec ular in 
a region of size a few hundred AU (jYoshida et al.l l2QQ6). 
Somewhat later, the gas becomes optically thick and the 
cooling radiation can no longer escape, inducing a final 
heating phas e that ends with the formation of a proto- 
star (Yoshida et al.l [2008). At this point a sink particle 
with an accretion radius of 100 solar radii is inserted, 
which is t he maximum photospheric s ize of the accreting 
protostar (jHosokawa fc Omukaill2009f ). 

Although the evolution of the gas is quite similar in 
all c ases, some scatter in th e thermal history exists 
(e.g., lO'Shea NormanI [2007t ). In simulations MH-2 
and MH-3, HD cooling becomes important and leads to 



a prolonged coolin g phase at d ensities nn ^ 10^ cm ^ 
(jRipamonti 2007; McGreer fc B rvan 2008). This is evi- 
dent from Figures 4 and 5, where we show the occupa- 
tion of the gas in density and temperature space. There 
is no correlation of the prolonged cooling with the bulk 
properties of the underlying DM halo (see Table 1), indi- 
cating that details of the virialization process must be re- 
sponsible for this effect. Note that the temperature floor 
required for HD cooling is around 150 K, which is only 
marginally attained in most minihalos. The activation 
of HD cooling is therefore very sensitive to the evolution 
of the gas during the initial free-fall phase. A second 
interesting effect is the shock-heating of some parcels of 
gas, which we attribute to the more violent heating in 
simulations MH-2 and MH-3 once three-body reactions 
set in. A peak similar to the one foun d at nn ^ 10^ cm~ ^ 
in simulation MH-2 was also found in lTurk et al.l (|2010f ). 
although it was much more pronounced. 

In Figure 6, we show the internal structure of the mini- 
halos just before the formation of the first protostar. The 
density, enclosed mass, and angular momentum profiles 
are relat ively similar and agree well with previo us investi- 
gations (jAbel et al.ll2002l : lYoshida et al.lbool ). The sec- 
ond density peak in simulation MH-2 indicates the nearly 
simultaneous colla pse of a seco nd gas cloud, which was 
also found in Turk et aTl ()2009f ). The temperature pro- 
files show significant scatter due to the previously men- 
tioned differences in the virialization of the individual 
minihalos. We find no convincing correlation of the ther- 
mal history of the gas with the velocity dispersion. The 
latter is always close to Mach numbers M 1, indicating 
transsonic turbulence. 

3.2. Fragmentation and Protostellar Mass Function 

Shortly after the formation of the first protostar, the 
residual angular momentum of the surrounding gas cloud 
leads to the formation of a circumstellar accretion disk 
(^Clark et al...2011b.) . Since the infall rate of new mate- 
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Fig. 14. — Resolution study of simulation MH-1 using standard sink particles. The top panels show the initial fragmentation of the disk 
in a box with 50 AU per side after 25 yr and for 8, 16 and 32 cells per Jeans length. The individual simulations are very similar during 
the initial formation and fragmentation of the disk, with the exception that in the middle panel two of the protostars have just merged. 
At later times the simulations begin to diverge quite strongly due to chaotic N-body interactions between the protostars. As shown in the 
bottom panels, the morphologies have become quite distinct after 200 yr, even though the total numbers of protostars are nearly equal. 



rial onto the disk is extremely high, angular momentum 
transport by viscous and gravitational torques is not suf- 
ficient to prevent the disk from fragmenting. As shown 
in Figure 7, the disk hence passes the gravitational in- 
stability limit, at which point it forms a number of sec- 
ondary protostars. After this initial burst, the gas be- 
comes gravitationally unstable multiple times and forms 
a small cluster of ~ 10 protostars after only 1000 yr (see 
Figure 8). 

The complex gravitational interactions between indi- 
vidual protostars and the surrounding gas clouds are il- 
lustrated in Figures 9 and 10, where we show the mass 
accretion histories of the protostars for two different 
sink-particle schemes. In our standard approach, sink 
particles merge if they become gravitationally bound 
and pass within 100 solar radii of each other, while 
in our 'adhesive' approach, only the latter criterion 
must be fulfilled. The second model is a crude rep- 
resentation of the dissipative effects of shocks in col- 
lisions between physically extended protostars. Both 
techniques have been used in (radiation) hydrodynamic 
simulations of star cluster formation in the present- 
day unive rse (Klessen & Burkert 2000; Bat e et al.l 120031: 
Krumholz et al. 2009; Federrath et al. 2010; Peters et al. 
20T0h . In the absence of a detailed understanding of the 
true protostellar radii and the influence of tidal forces 
during close encounters, we hope that they bracket phys- 
ical reality. 

Interestingly, the two methods predict nearly identi- 



cal total protostellar masses, but yield different detailed 
mass accretion histories. For standard sink particles, 
protostars generally survive close encounters, whereas for 
adhesive sink particles these interactions often result in 
a merger. For this reason the number of protostars af- 
ter 1000 yr is higher in our standard formulation of sink 
particles than in our adhesive formulation. This is evi- 
dent from Figure 11, where we show the final protostel- 
lar mass functions of all five minihalos. Notably, a small 
group of protostars with a range of masses is formed in 
both cases. The cumulative protostellar mass functions 
obtained from the five independent realizations are rel- 
atively flat, extending from ~ 0.1 to nearly 10 Mq 
(bottom right panel of Figure 11). This corresponds to 
a top-heavy mass function, in the sense that most of the 
protostellar mass is locked up in massive protostars. 

An interesting trend obtained in the case of stan- 
dard sink particles is that some low-mass protostars are 
ejected from the central gas cloud by dynamical interac- 
tions with more massive protostars, while the latter tend 
to remain at the center of the minihalo and continue to 
accrete aggressively. This is also visible in Figure 12, 
where we show the radial velocity and ratio of the ra- 
dial velocity to the escape velocity of all protostars. A 
number of protostars move away from the center of mass 
with speeds exceeding the escape velocity, such that they 
stop accreting and might enter the main sequence as low- 
mass Pop III stars. This still occurs in our formulation 
of adhesive sink particles, but less frequently. Further- 
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Fig. 15. — Same as Figure 9, but for the resolution study of simulation MH-1. The three runs show some differences in the mass accretion 
histories of individual protostars due to the chaotic nature of N-body interactions. However, their numbers and total mass, as well as their 
cumulative and average accretion rates agree well with each other. The employed resolution therefore does not qualitatively affect the 
fragmentation of the gas. 



more, in this case no protostars reside within the central 
100 AU, since momentum conservation acts to move 
merged protostars to larger annuli. 

3.3. Radiation Feedback 

Studies of present-day star formation have shown that 
the radiation emitted by young protostars heats the sur- 
rounding envelope and reduces t he degree of fragmen- 
tation near the c entral protostar (iKrumholz et al. 2009^ 
iPeters et al.ll2010l) . However, the outcome is less clear in 
the present study, where gas temperatures are a factor of 
10 — 100 times higher and different cooling mechanisms 
are in play. To investigate the importance of radiation 
in detail, we compare simulations with and without the 
heating from accreting Pop III stars in one of the five 
minihalos. 

During the early stages of protostellar evolution. 
Pop III stars grow rap i dly in size from :^ 10 to 100 R© 
(|Hosokawa fc Omukail 120091 ). The temperature of the 
photosphere remains roughly constant at 6000 K, and so 
most of the energy is radiated in the optical and infrared. 
Once the protostar has accreted of order 10 M©, the star 
undergoes Kelvin-Helmholtz contraction and heats up 
substantially, leading to the emission of ionizing radia- 
tion. Since we are well within the initial swelling phase, 
we can utilize a simple optically thin treatment of the ra- 



diation to determine the maximum extent to which the 
gas will be affected. 

For this purpose we treat each sink particle formed 
in our simulation as a separate protostar, and account 
for the energy released by accretion onto the surfaces of 
these protostars. To model the effects of this accretion 
luminosity, we first compute the bolometric accretion lu- 
minosity for each protostar 



(9) 



where is the accretion rate onto the protostar, 
is the protostellar mass and is the protostellar ra- 
dius. We relate the protostellar radius to the cur- 
rent protostellar mass and accretion rate using a rela- 
tionship derived for adiabatically accreting, metal-free 
protostars eni bedded in a spherically symmetric infiow 
(jStahler et anil986l) : 



i?.=26Re(^ 



0.27 



0.41 



10-3 Mo yr-i 



(10) 



More recent calculations for the case of rotating infall 
find an even large r value f or R^ during the adiabatic 
accretion phase (^Tan fc McKee„.2004) . and so using the 
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Fig. 16. — The final mass function obtained in simulation MH-1 using accretion radii decreased by factors of two and four from the 
fiducial value race = 100 R©. The collapse of the gas is thus followed to higher densities and smaller separations from the protostar. Any 
systematic inaccuracies stemming from the artificially reduced density and pressure around sink particles should therefore induce a trend. 
As is evident from the figure, such a trend exists: the amount of fragmentation increases, while the typical mass of a fragment decreases. 
This shows that our results should be considered a lower limit on the degree of fragmentation in minihalos. 



above expression gives us a conservative upper limit on 
the bolometric accretion luminosity during this phase. 
Combining the above equations, we find that 



1200 L 







0.73 



0.59 



10-3 Mq yr- 



(11) 

With this expression, we determine the heating rate of 
the gas surrounding the protostar from 



47rr^ 



(12) 



where p is the mass density, r is the distance to the pro- 
tostar, and is the Planck mean opacity of the gas. We 
calculate this mean opacity by interpolation, using tab- 
ulated valu es that account for bot h line and continuum 
absorption (jMaver fc Duschill2005f ). and that include the 
influence of the electrons provided by ionized lithium. 
This expression assumes that the gas is optically thin to 
the radiation from the accreting protostar. Making this 
assumption allows us to avoid the extremely high compu- 
tational cost that would be associated with an accurate 
treatment of the transfer of the protostellar radiation, 
and also gives us a conservative upper limit on the effec- 
tiveness of protostellar feedback. Finally, we assume that 
each protostar accretes at a rate of 0.1 M© yr"^, which is 
well above the total accretion rate of the whole ensem- 
ble of protostars, maximizing the effect of our radiative 
feedback model. 

In Figure 13, we compare runs with and without ra- 
diation feedback for simulation MH-1. We flnd that the 
additional heating results in a smoother and more ex- 
tended disk, in which the fragments form at slightly 
larger distances from the center (see also iClark etldl 
[2011bh . However, even the unrealistically strong radi- 
ation field assumed here does not prevent fragmenta- 
tion within the disk. The primary reason for this is the 
strong cooling of the gas by molecular hydrogen lines at 
temperatures around 1000 K and above, which largely 
offset s the heating cause d by the accretion luminosity 
(e.g., lOmukai et al.l 12010^. The neglect of radiation in 
our simulations should therefore not qualitatively affect 



our conclusions. We note that a more detailed investi- 
gation of radiation feedback for the minihalos used her e 
is presented in a companion paper ([Smith et al.l 120111 ). 
Finally, it is important to point out that the heating 
may be much more severe once the first protostar has 
accrete d of order 10 M,^ and begins to emit ioniz ing ra- 
diation (|^^EMcKei[200l |McKeeX^3 How- 
ever, this will not occur until after the period of time 
simulated here. 

3.4. Resolution Study 

We here study how our results depend on the reso- 
lution employed. For this purpose we rerun simulation 
MH-1 with a modified refinement criterion. We linearly 
interpolate the number of cells per Jeans length between 
two threshold densities, using Jinit = 1/128 at a density 
of nn = lO^cm-3 and Jfinai at nn = 10^ ''cm^^. This 
allows a better control of the number of resolution ele- 
ments at a given density, and thus on the scale where the 
sink particles are created. We run three simulations with 
inverse Jeans numbers of Jfinai = 1/8,1/16 and 1/32, 
corresponding to factors of 8 increase in the number of 
particles at nn = 10^^ cm"^. In analogy to our main 
simulations, sink particles are created once Jcrit = Jfinah 
with the exception that we here use an initial sink par- 
ticle radius equal to the accretion radius of 100 R©. 

The mass resolution of these simulations is significantly 
higher than in our original simulations, which greatly in- 
creases their computational cost. The mass in the cells at 
the highest densities has dropped by nearly two orders 
of magnitude to lO^^M©, so that we have only in- 
vestigated the first 200 yr instead of the original 1000 yr. 
However, this period of time is sufficient to show that 
the overall fragmentation process is not affected by the 
resolution employed. In Figure 14, we show the central 
50 AU of all three simulations after 25 yr. The morpholo- 
gies of the disk and the positions of the protostars agree 
well with each other, with the exception that in the case 
with Jfinai = 1/16 one of the protostars has just merged. 
Somewhat later, N-body interactions lead to an increas- 
ingly chaotic evolution of the system, so that small differ- 
ences in the initial conditions result in a relatively large 
dispersion in the evolutionary paths of individual proto- 
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stars. This is evident from the bottom panels, where we 
show the protostellar cluster in a box with 200 AU on 
a side and after 200 yr. The appearance of the cluster 
has become quite distinct in each case. However, Fig- 
ure 15 shows that the numbers of protostars and their 
total mass are much more insensitive to the initial con- 
ditions. The employed resolution can thus change the 
details of the fragmentation due to a 'butterfly effect', 
but not qualitatively affect our results. 

3.5. Variation of Accretion Radius 

The removal of mesh-generating points below the ac- 
cretion radius artificially reduces the density and pres- 
sure around sink particles. We investigate how this af- 
fects our results by systematically decreasing the accre- 
tion radius from 100 R© to 50 R© and 25 R©. The gas 
therefore collapses to progressively higher densities be- 
fore being accreted onto the sink particle, such that any 
residual effects should display a trend in the total num- 
ber of fragments formed, as well as on the masses of the 
fragments. As is evident from Figure 16, such a trend 
exists: A reduction of the accretion radius results in an 
increase in the number of fragments and a decrease of 
the typical fragment mass. Existing sink particles ap- 
parently accrete less aggressively, such that more mass 
remains in the disk, which in turn becomes even more 
susceptible to fragmentation. This trend indicates that 
our results should be considered a lower limit on the de- 
gree of fragmentation in minihalos. 

3.6. Protostellar Merging 

Another limitation of the sink-particle algorithm em- 
ployed here is that the gasdynamical friction between 
physically extended protostars is not modeled. Mergers 
between protostars may therefore be more common than 
for the purely gravitationally interacting sink particles 
employed here. To understand the potential importance 
of gasdynamical friction, we use Equation 10 to compare 
the minimum separation of each protostar to any other 
protostar with the sum of the protostellar radii during 
their closest encounter. Figure 17 shows that the major- 
ity of protostars experience at least one close encounter 
with another protostar and, depending on the eccentric- 
ity and infall velocity of the respective orbit, will merge. 
Only a few low-mass protostars remain physically sepa- 
rated for the entire duration of the simulation, while ac- 
quiring substantial radial velocities (see Figure 12). As 
mentioned already in Section 3.2, these protostars stop 
accreting and may therefore enter the main sequence as 
low-mass Pop III stars. However, this possibility depends 
sensitively on our choice of the accretion radius, which 
might be too small considering that real protostars are 
likely surrounded by an unresolved accretion disk. As- 
suming a somewhat larger accretion radius would further 
reduce the number of low-mass protostars that could 
survive. We note that an accurate account of merg- 
ing between protostars would require simulations that 
self-consistently capture the interaction of the protostars 
with the surrounding gas cloud. These are not currently 
feasible. 

4. SUMMARY AND CONCLUSIONS 

We have used the moving mesh code AREPO to fol- 
low the runaway collapse of the gas in five statistically 
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Fig. 17. — A comparison of the minimum separation of each pro- 
tostar to any other protostar with the sum of the protostellar radii 
during their closest encounter. Black dots, crosses and stars de- 
note protostars with masses below 1 Mq, between 1 Mq and 3 M©, 
and above SM©. Nearly all protostars experience at least one en- 
counter with another protostar where the separation drops below 
their total physical size. Depending on the eccentricity and infall 
velocity of the respective orbits, they will merge. As delineated 
by the dotted line, only a few low-mass protostars remain physi- 
cally separated for the entire duration of the simulation. We note 
that this number depends sensitively on our choice of accretion 
radius, which might be too small considering that real protostars 
are surrounded by an unresolved accretion disk. Assuming a some- 
what larger accretion radius would further reduce the number of 
low-mass protostars that could survive. 



independent minihalos from cosmological to protostellar 
scales - over more than twenty orders of magnitude in 
density. We have captured the subsequent evolution of 
the newborn protostellar cloud for more than 100 dynam- 
ical times with a sink-particle algorithm that resolves the 
accretion of the gas down to scales of 100 R© or less, com- 
parable to the maximum photospheric size of Pop III pro- 
tostars. As proposed bv IClark et al.l (|2011b[ ). in all five 
minihalos a circumstellar disk forms that fragments vig- 
orously into a small cluster of protostars with a range of 
masses. The gas becomes gravitationally unstable multi- 
ple times and continues to form protostars for the entire 
duration of the simulation. After only 1000 yr of accre- 
tion, of order 10 protostars per minihalo have formed 
and display a relatively flat protostellar mass function 
ranging from ~ 0.1 to nearly 10 M©. 

Although the sink-particle technique employed in this 
study allows a continuation of the calculations to much 
later times than was possible in previ ous ab initio sim- 
ulations of primord ial star formation ( !Abel et all 120021 : 
lYoshida et al.l [20081 ), a number of uncertainties remain. 
One problem is the artificially reduced density and pres- 
sure created by the removal of mesh-generating points 
around sink particles. However, we have found that us- 
ing systematically smaller accretion radii results in more 
fragmentation and a decrease of the typical fragment 
mass, such that a correct treatment of the boundary re- 
gion between sink particles and the gas should strengthen 
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our conclusions. A more important caveat is the inability 
of the sink-particle algorithm to model the gasdynamical 
friction between physically extended protostars during 
close encounters. In an attempt to maximize this effect, 
we have implemented adhesive sink particles in addition 
to a more standard formulation of sink particles. This 
leads to an increased merger rate, so that fewer proto- 
stars survive. However, even this extreme assumption 
does not prevent the formation of a relatively rich clus- 
ter of protostars. A final caveat is that we do not self- 
consistently model the interaction of the protostars with 
the surrounding gas. It is unclear how important the re- 
sulting neglect of gasdynamical friction and torques is, 
since simulations that model the protostellar surface as 
well as the parent gas cloud are not yet feasible. However, 
it appears unlikely that this caveat will qualitatively af- 
fect our conclusions, since fragmentation typically occurs 
on scales larger than the minimum resolution length. 

Modulo the uncertainties mentioned above, the sim- 
ulations presented here portray a very different picture 
of primordial star formation than is commonly assumed. 
Instead of forming a single object, the gas in minihalos 
fragments vigorously into a number of protostars with a 
range of masses. It is an open question as to how this 
early mass function will be mapped into the final mass 
function of Pop HI, after accretion, fragmentation and 
merging have finally stopped. However, it is interest- 
ing to speculate how this nonlinear mapping will play 
out. If a fiat, broad mass function persists, a number of 
lower-mass Pop HI stars will have formed which could 
survive to the present day if their mass remains below 
~ O.SMq. Although this possibility is speculative be- 
cause of the above uncertainties and the fact that we 
follow the protostellar accretion only for the first 1000 
out of 10^ or 10^ yr, it is worth looking for such Pop HI 
fossils in ongoing and planne d larg:e surveys of metal- 
poor stars in the Milky Way (iBeers fc Christliebl 12005). 
Specifically, such surveys should be focused on the Galac- 
tic bulge, where the Pop HI survivors should preferen- 
tially reside due to the biasing of the min ihalo forma- 
tion sites (|Diemand et~al] 120051 : iGao et al.l 1201 0). The 
planned Apache Point Observatory Galactic Evolution 
Experiment (APOGEE) with its near-IR capa bility may 
be well suited for this search (' Majewski et al.ll2010>) . Be- 
cause of interstellar pollution, even true Pop HI fossils 
would appear as extreme Pop II stars, but upper limits 
would be significantly low er than the currently probed 
values (jPrebel et al°]l2009f ). 



Furthermore, the presence and mutual competition of 
multiple accretors in a given minihalo will act to limit 
the growth of the most massive objects. Pop HI stars 
are therefore less likely to reach masses in excess of 
^ 140 M0, the threshold for triggering extremely ener- 
getic p air-instability supernova e (PISNe) during stellar 
death (jHeger WooslevI l2002f ). A reduced PISN rate 
is more easily compatible with the absence of their dis- 
tinct nucleosynthetic signatures in any o f the extremely 
metal-poor halo stars observed so far (jlwamoto et al.l 
l2005l ). Pop HI stars could still have given rise to nu- 
merous extremely luminous supernova explosions, if they 
had masses of a few 10 M© and if th ey were rapid rota - 
tors, as suggested by recent studies (jStacv et al.ll201lh . 
They would then have exploded as core-collapse hy- 
pernova e, with explosion energi es that are similar to 
PISNe (jUmeda Nom oto' '20021. Finally, our results 
may challenge models of so-called 'dark stars', which are 
Pop HI stars powered by DM self-a nnihilation heating 
(|Freese et al.l 120081 : llocco et al.ll2008f ). These models in- 
voke an increased DM interaction rate at the center of 
the Pop HI star which itself lies at rest at the center of 
its minihalo. On top of recent claims that the initial col- 
lapse is not substantia lly delayed by DM annihilations 
(jRipamonti et al.]l2010l ). the complex dynamics of a pro- 
tostellar cluster at the center of the minihalo may further 
preclude efficient DM capture and heating. 
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